Multiple-relaxation-time lattice Boltzmann model for 

compressible fluids 

Feng Chen^, Aiguo Xu^*, Guangcai Zhang^, Yingjun Li^ 
1, China University of Mining and Technology (Beijing), Beijing 100083 
2, National Key Laboratory of Computational Physics, 
Institute of Applied Physics and Computational Mathematics, 
P. O. Box 8009-26, Beijing 100088, P. R. China 
(Dated: April 20, 2010) 

Abstract 

We present an energy-conserving multiple-relaxation-time finite difference lattice Boltzmann 
model for compressible flows. This model is based on a 16-discrete-velocity model. The colli- 
sion step is first calculated in the moment space and then mapped back to the velocity space. The 
moment space and corresponding transformation matrix are constructed according to the group 
representation theory. Equilibria of the nonconserved moments are chosen according to the need of 
recovering compressible Navier-Stokes equations through the Chapman-Enskog expansion. Numer- 
ical experiments showed that compressible flows with strong shocks can be well simulated by the 
present model. The used benchmark tests include (i) shock tubes, such as the Sod, Lax, Colella 
explosion wave, collision of two strong shocks and a new shock tube with high mach number, 
(ii) regular and Mach shock reflections, (iii) shock wave reaction on cylindrical bubble problems, 
and (iv)Couette flow. The new model works for both low and high speeds compressible flows. 
It contains more physical information and has better numerical stability and accuracy than its 
single-relaxation-time version. 
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I. INTRODUCTION 



The Lattice Boltzmann (LB) method is an innovative numerical scheme originated from 
Lattice Gas Automata (LGA) l| and aim to simulate various hydrodynamics |2|. The LB 
method was introduced to overcome some serious deficiencies of LGA, such as intrinsic 
noise, limited values of transport coefficients, non-Galilean invariance, and implementa- 
tion difficulty in three dimensions. In the past two decades, most of the LB models are 
based on the famous Bhatnagar-Gross-Krook (BGK) approximation [sl] where a Single Re- 
laxation Time (SRT) is used. Due to its validity and simplicity, the SRT LB method has 
been widely used to simulate various fluid flow problems, such as the multiphase flowjj- 
magnetohydrodynamics 12Nl8l|. flows through porous media 19N21| and thermal fluid 



dynamics 



22- 



2^, etc. 



However, the extreme simplicity of the SRT leads also to some constraints for the SRT 
LB model. For example, the simulation will be unstable when the relaxation time r is 
close to 0.5, the model works only for low Mach number flows. One possible remedy is to 
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26| . In real fluid the equilibrating 



use the Multiple Relaxation Time (MRT) method 
rates of mass, momentum, energy, etc. are generally different. This difference can be 
manifested by the non-unique adjustable parameters in the MRT LB model. In contrast 
to the SRT model, the MRT version has much more adjustable parameters and degrees 
of freedom. The relaxation rates of various processes owing to particle collisions may be 
adjusted independently. The main strategy of the MRT LB scheme is that the collision step 
is flrst calculated in the moment space and then mapped back to the velocity space. The 
advection ste p is still computed in the velocity space. In many cases, it has been shown by 
Luo, et al 
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that the MRT LB model has better numerical stability. 



Recently, the MRT LB method has attracted considerable interest and much progress 



nas 



been achieved. For example, MRT models for viscoelastic fluids 
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-|31|. multiphase flows 
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33| . flow with free surfaces 34l| . etc. were developed; optimal boundary condition for MRT 



LB was composed |35|. To simulate system with temperature fleld, Luo, et al.[28| suggested 
a hybrid thermal MRT LB model. These models work only for nearly incompressible fluids 
with very low Mach number. 



LB community has long been attempting to construct models for compressible fluids 
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39 |. Alexander and Chen et al. j36| constructed a model where the sound speed is adjustable 
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so that the Mach number can be enhanced. Li, et ah 37|] gave a model by reforming the 



velocity space. Sun, et al. [SSj formulated adaptive LB models where the particle velocities are 
determined by the mean velocity and internal energy. Yan, et al. [39] proposed three-speed- 
three-energy-level models. Besides the standard LB mentioned above, some researchers have 
also tried to develop Finite Difference (FD) LB for compressible fluids 4Ch42| . but in the real 
simulations the accessible Mach number is still not large. The model introduced by Kataoka 



4l| uses only sixteen discrete velocities and hence has a high computational 



and Tsutahara 
efficiency. 

The low-Mach number constraint is generally related to a numerical stability problem. 
The latter has been partly addressed by a number of techniques, such as the entropic 



method 



431, 



44j |. the fix-up scheme 43|, |45|, Flux-limiters[46!] and dissipation 47|, |48|] tech- 



niques. In existing SRT models, it seems that the most effective solution to overcome the 
low Mach number constraint is to introduce artificial viscosity. But with the artificial vis- 
cosity, some fundamental kinetics are not very clear. In many cases, the MRT formulation 
has been shown to offer improved numerical stability, and provide additional physics. In 
this paper we present an energy-conserving multiple relaxation time finite difference lattice 
Boltzmann model for compressible flows with high Mach number. This model is based on the 
one proposed by Kataoka and Tsutahara 4l|. The moment space and transformation matrix 
are constructed according to the group representation theory. Equilibria of the nonconserved 
moments in the moment space are chosen when recovering compressible Navier-Stokes (NS) 
equations through the Chapman- Enskog (CE) expansion. 

This paper is organized as follows. In Sect. II a brief review to the MRT LB model is 
presented. In Sect. Ill the new model is constructed. The von Neumann stability analysis is 
given in Sect. IV. Section V shows the numerical tests and some simulation results. Section 
VI provides a summary and concludes the paper. 



II. BRIEF REVIEW OF THE MRT LB MODEL 



The evolution of the distribution function f] for the particle velocity Vi is governed by 
the following equation: 



dt 
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where fi (/f^) is the particle (equihbrium) distribution function, Vi represents a group of 
particle velocities, subscript i indicates the particle's direction, i = l,...,N,Nis the number 
of discrete velocities, the subscript a indicates x or y component, S is the collision matrix. 
The equation reduces to the usual lattice BGK equation if all the relaxation parameters are 
set to be a single relaxation time r, namely S = ^I, where I is the identity matrix. 

The discrete (equilibrium) distribution function fi {f^'^) in Eq. ([1]) can be listed with the 
following matrixes: 

f=(/i,/2,---,/7vr, (2a) 
= (/r, , ■ ■ ■ , fN'f , (2b) 

where T is the transpose operator. 

Given a set of discrete velocities f j, and corresponding distribution functions /j, we can get 
a velocity space 5*^, spanned by discrete velocities f j, and a moment space S'*'^, spanned by 
moments of particle distribution function /j, where i = 1, - ■ ■ ,N . Similarly, we also express 
the moments of distribution function with the column matrix: f = ^/i, /2, ■ ' ' ; /a? j ? where 
fi = rriijfj, rriij is an element of the matrix M and is a polynomial of discrete velocities. 
Obviously, the moments are simply linear combination of distribution functions fi, and the 
mapping between moment space and velocity space is defined by the linear transformation 
M, i.e., f = Mf , f = M^^'^f , where M = (mi, m2, ■ ■ ■ , m^v)^ , rrii = (m^i, mj2, ■ ■ ■ , ran^). 

The LB simulation consists of two steps: the collision step and the advection one. In the 
MRT LB method, the advection step is computed in the velocity space. The collision step is 
first calculated in the moment space and then mapped to the velocity space. So, the MRT 
LB equation can be described as: 

^+^.a^ = -M-iS,„(A-/f), (3) 

where S = MSM^^ = diag(si, S2, ■ ■ ■ ,sn) is the diagonal relaxation matrix, f^'^ is the 
equilibrium value of the moment fi. The moments can be divided into two groups. The first 
group consists of the moments locally conserved in the collision process, i.e. fi = f^^. The 
second group consists of the moments not conserved, i.e. fi ^ f^''. The equilibrium f^'' is a 
function of conserved moments. 
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FIG. 1: Distribution of ^^ia for the discrete velocity model. 



III. ENERGY-CONSERVING MRT LB MODEL 



We use the two-dimensional discrete velocity model by Kataoka and Tsutahara 4l| (see 
Fig. 1). It can be expressed as: 

' eye : (±1, 0) , for 1 < i < 4, 
eye : (±6, 0) , for 5 < i < 8, 
v^(±l,±l), for9<i<12, 
^ ^(±1,±1), for 13<z<16, 



(4) 



where eye indicates the cyclic permutation. 



A. Construction of transformation matrix M 

The transformation matrix M is constructed according to the irreducible representations 
of SO (2) group: 

1, 

cos 6, sin 0, 

sin^ 6 + cos^ ^, cos 2^, sin 26', 

cos 6'(sin^ 6 + cos^ 6*), sin 6'(sin^ 6 + cos^ 6*), cos 36', sin 36', 

(sin^ e + cos^ ef, cos 2e{si\i^ 9 + cos^ 9), sin 2^(sin2 9 + cos^ 9), cos 40, sin 40, 
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Let Vix and Viy play the roles of cos^ and sin^, respectively. Then we define 



mii = 1, 


(5a) 


m2i = Via:, 


(5b) 




(5c) 




(5d) 




(5e) 




(5f) 




(5g) 




(5h) 




(5i) 


mioi^Viy{3vl-v^y), 


(5j) 


^iii = (Vix + ViyY/A, 


(5k) 


mi2i = vf^ - Qv^v^y + vfy, 


(51) 


mm = {vl + vl){vl - wj), 


(5m) 


^i4i = {vfx + Viy)VixViy, 


(5n) 


mi5i -= Vix{vl - 3vl){vl + v^), 


(5o) 


miQi = Viy{2>vl^ - vl){vl^ + vl). 


(5p) 



where i = 1, ■ ■ ■ ,16. 

For two-dimensional compressible models, we have four conserved moments, density p, 
momentums jx, jy, and energy e. They are denoted by /i, /2, fs and f^, respectively. 
Specifically f^^p^ h = jx = E/i"^2i, /s = jy = Y^fimsi, /4 = e = Y.fi'^n- 

To be consistent with the idiomatic expression of energy, in the definitions of m^i, niji and 
mgi, a coefficient 1/2 is used. Similarly, a coefficient 1/4 is used in the definition of mm. 
Thus, the transformation matrix M can be expressed as follows: 

M = (mi, 7712, • • • ,mie)'^, 
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where 



1712 



mi = (l, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1,1, 1,1), 
(1, 0, -1, 0, 6, 0, -6, 0, V2, -V2, -V2, ^2, -^), 

(0, 1, 0, -1, 0, 6, 0, -6, V2, V2, -V2, -V2, --^), 

,1 1 1 1 9 9 9 9, 

m4 = (-,-,-,-,18,18,18,18,2,2,2,2,-,-,-,-), 
* ^2222 '''''222 2^ 

ms = (1, -1, 1, -1, 36, -36, 36, -36, 0, 0, 0, 0, 0, 0, 0, 0), 

me = (0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 2, -2, ^, ^, -^), 

1 1 27 27 27 27 

= (_ 0. 0. 108. 0, -108, 0, 2V2. -2V2. -2^2. 2^2. ^. 

1 1 27 27 27 27 

ms = (0, -, 0, --, 0, 108, 0, -108, 2^2, 2v^, -2v^, -2v^, -^), 

27 27 27 27 

mg = (1, 0, -1, 0, 216, 0, -216, 0, -4^2, 4^2, 4^2, -4^2, -, — , — , ^ 



27 27 27 27 

mio = (0, -1, 0, 1, 0, -216, 0, 216, 4^2, 4v^, -4v^, -4v^, — , — , - — , -— ), 

= ^i' i' i' i' ^^^' ^^^' ^^^' ^^^' ^' ^' ^' ^' T' f ' T' T^' 
mi2 = (1, 1, 1, 1, 1296, 1296, 1296, 1296, -16, -16, -16, -16, -81, -81, -81, -81), 

mi3 = (1, -1, 1, -1, 1296, -1296, 1296, -1296, 0, 0, 0, 0, 0, 0, 0, 0), 

mi4 = (0, 0, 0, 0, 0, 0, 0, 0, 8, -8, 8, -8, y , -y , y , -y ), 

^ ^ ^ ^ 243 243 243 243, 

mi5 = 1, 0, -1, 0, 7776, 0, -7776, 0, -16 A 16^2, 16^2, -16^2, ^, ^ , 

W2 V2 V2 V2 

r r r- ^ 243 243 243 243, 

mi6 = 0, -1, 0, 1, 0, -7776, 0, 7776, 16v^, 16^2, -16^2, -16^2, ^, ^ . 

72 72 ^2 W2 

B. Determination of /^"^ 



We perform the Chap man- Enskog expansion [321. |49|, |50|] on the two sides of Eq.([T]). We 
define 

h = ff'^ff'^f?. (6a) 
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where is the zeroth order, f-^\ ^and ^ are the first order, /-^^^ and ^ are the second 
order terms of the Knudsen number e. Equating the coefficients of the zeroth, the first, and 
the second order terms in e gives 

/f = ff, (7a) 

They can be converted into moment space to obtain: 

ff' = ff. (8a) 

where E^, = M{vial)M-'^. 

The equihbria of the moments in the moment space can be defined as : f^'' = 
{p, jx, jy, f^'^ , /q"^ , ■ ■ ■ , /le)^- The first order deviations from equilibria are defined as : 
f^^^ = {0,0,0,0, f^^\ fQ^\ ■■ ■ j/ie'')"^. In the same way, the second order deviations are 
f(2) = (0, 0, 0, 0, ff\ fP, ■■■ , fiff. From Eq.(l8b]) we obtain 

*+|^ + §i = o. (9a) 

oti oxi oyi 

|^r+^(/s+i/g) + i|^/s^-«#> (9g) 
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/,'!', (9k) 



|/J'4^(/3 + /S)-|^/S^Wi'> (9i) 

+ + = -snfil (9m) 

^/d + ^^i-Qjy + \ft! + §fl'-^^^^^^ = -SufiW (9n) 

|^AH^(75e-54p-18/f +18/rf+|/r|+y/r|) + ^(36/f ^ = s^J^l (9o) 



From Eq. llHcI) we obtain 



1^0. (10a) 



Adding Eq.(lO) and the first four formulas of Eq.(9) leads to the following equations, 

^ + ^ + ^ = dla) 
dt dx dy 



djx 
dt 








1 d ^(1) .9 ^(1) 
2dx^^ dy^'^ ' 


(lib) 


dt 


+ ±f'- + ^1 

dx ^ dy 






^ f (1) I 1 ^ f (1) 


(11c) 




^ + 9 + 
dt dx ^ 


dy^' 


d 
dx 


P(i) 5 


(lid) 



we choose 



f^' = {£-£)/P, (12a) 
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/r = j.jy/p, (12b) 

f^'^ = (e + P)jjp, (12c) 

/f = (e + P)jy/p, (12d) 

/r = (i^ - ^31)3./P\ (12e) 

m = 2e'/p-{3l+jlflAp\ (12g) 

Al = (6pe - 2jl - 2jl){jl - jl)/p\ (12h) 

fit - (6pe - 2jl - 2jl)j,jjp\ (12i) 

The definitions of /^f, /fl, fH have no effect on macroscopic equations, so we can choose 
/i2 = fit = /i6 ~ 0- l^ ^^i^ "^^y ^be recovered NS equations are as follows: 

^ + ^ + ^ = 0, (13a) 
dt dx dy 

. d ,^ . d dP d du^ duy d duy du^ 

^iy , "9 (9 ..2, X dP d duy du^ d du^ duy 

| + ^[(e + P),y,] + |[(e + P),yp] 

(9 dT duy dux duy du^,. 



9x dx ^ dx ^ dx ^ dy ^ dy 

where //^ = pRT/s^, py — pRT/sQ, Ai = pRT/s-j, A2 = pRT/ss- When ps — Pv — P, 
Ai = A2 = A, the above NS equations reduce to 



■'a 



^ ^ '9 (iaj/g/p) ^ _ ^ + ^ _ )1 (14b) 

dt dxp dxa dxp dxp dxa dx^ 

de d . , T d r^,^dT ,dua dug duy ^ , , 
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IV. VON NEUMANN STABILITY ANALYSIS 



In this section we perform the von Neumann stabihty analysis on the new MRT LB 
model. In the stability analysis, we write the solution of FD LB equation in the form of 
Fourier series. If all the eigenvalues of the coefficient matrix are less than 1, the algorithm 
is stable. 

The distribution function fi{xa,t) is split into two parts: 

Mx^,t) = f^ + AMx^,t), (15) 

where is the global equilibrium distribution function. It is a constant and does not vary 
with time or space, depends only on the average density, velocity and temperature. Similarly, 
the distribution function fi{xa,t) is split into two parts: 

U^a,t) = f) + AUx^,t), (16) 

where /° is a constant. Putting the Eq. ffTSl) and Eq. fflGl) into Eq.([3]) gives 

AMx^,t + At) = Af,ix^,t) - Atv,^^ - AtMT^'SikiAfk - A/f ). (17) 

OXa 

The solution can be written as 

Afi{xa,t) = F^exp{ikaXa), (18) 

where F* is an amplitude of sine wave at lattice point Xa and time t, is the wave number. 
From the above two equations we can obtain _p^*+^* = djFj. Coefficient matrix Gij describes 
the growth rate of amplitude F/ in each time step At. If u denotes the eigenvalue of 
coefficient matrix Gij, the von Neumann stability condition is max|a;| < 1. Coefficient 
matrix G^j can be expressed as follows, 

+ e-''^"^-)5., - AtM,^S,.(|^ - ^), (19) 



where 



fk — ^kpfp, 

df, dp dfj dT dfj dUa df, ■ ^ ^ 
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FIG. 2: Stability comparison for the new MRT model and its SRT version. 

Coefficient matrix Gij contains a large number of matrix elements (as many as 16 x 16), 
and every matrix element correlates with the macroscopic quantities and model parameters, 
so analytic analysis is very difficult. We conduct a quantitative analysis using Mathematica 
software. 

In Fig. 2 we show an example of stability comparison for the new MRT model and its 
SRT version. The abscissa is for kdx^ and the vertical axis is for Iwlma^: which is the largest 
eigenvalue of coefficient matrix Gij. The macroscopic values in stability analysis arc chosen 
as follows: {p,ui,U2,T) — (2.0,10.0,0.0,2.0), other common parameters are: dx — dy — 
2 X 10~^, dt — 10~^. The relaxation time in SRT is r = 10~^, while the coUision parameters 
in MRT are S5 = 6500, 57 = sg = 9 x 10^, sg = 8 x 10^ S13 = 7 x 10^, S14 = 8 x 10^ 
Si5 = 2.5 X 10^, the others are 10^. In this case, the MRT scheme is stable, while the SRT 
version is not. It is clear that, by choosing appropriate collision parameters, the stability of 
MRT can be much better than the SRT. 
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V. NUMERICAL SIMULATIONS 



In this section we study the following problems using the new MRT LB model: One- 
dimensional Riemann problems, shock reflections, shock wave reaction on cylindrical bubble, 
and Couette flow. 

A. One-dimensional Riemann problems 

Here, we study several typical one-dimensional Riemann problems, including the Sod, 
Lax shock tube, Colella explosion wave, collision of two strong shocks and a new shock tube 
with high Mach number. In the x direction, fi — f^*^ is set on the boundary nodes before 
the disturbance reaches the two ends. In the y direction, the periodic boundary condition 
is adopted. In the following part, subscripts "L" and "R" indicate the macroscopic variables 
at the left and right sides of discontinuity. 

(a) Sod shock tube problem 

The initial condition is 

(p,Mi,M2,T)U = (1.0,0.0,0.0,1.0), x<Q. ^^^^ 
{p,Ui,U2,T)\r = (0.125,0.0,0.0,0.8), x > 0. 

Figure 3 shows a comparison of the MRT LB simulation results and exact solutions for the 
density, pressure, velocity and temperature of the Sod shock tube problem at time t — 0.18. 
Here, in the collision matrix = sq = 5 x 10^, sj = sg = 10^, Sn = 2500, and other collision 
parameters arc 10^. The red circles correspond to simulation results with the grid size dx = 
dy — 0.002 and time step dt = 2x 10~^ (case 1), the green triangles correspond to simulation 
results with dx — dy — 0.001, dt — 10~^ (case 2), and solid hues represent the exact solutions. 
The relative errors of the density, pressure, velocity and temperature for case 1 are 0.234%, 
0.182%, 3.32% and 0.327%, respectively. The relative errors of the density, pressure, velocity 
and temperature for case 2 are 0.225%, 0.171%, 3.16% and 0.322%, respectively. Here, the 
relative error is defined as = V . k, , — "Jr , I / f kr , I , where , denotes 
the variables at the node of {xi,yj) obtained from the numerical simulation, and <Jjjg^„ is 
the exact solution at the same node. The simulation results successfully capture the main 
structure of the waves. 

b) Lax shock tube problem 

13 




FIG. 3: Comparison of the MRT LB simulation results and exact solutions for Sod shock tube at 
time t = 0.18. 

The initial condition of the problem is: 

(p, Ml, M2, r)U = (0.445,0.698,0.0,7.928), x < 0. ^^^^ 
{p,ui,U2,T)\r = (0.5,0.0,0.0,1.142), x > 0. 

Figure 4 shows the MRT LB numerical results and exact solutions for the Lax shock tube 
problem at time t = 0.2. The red squares, green circles and blue triangles correspond to 
simulation results with different grid sizes and time steps: dx = dy = 0.004, dt = 4 x 10^^ 
(case 1), dx = dy = 0.002, dt = 2 x IQ-^ (case 2), dx = dy = 0.001, dt = 10"^ (case 
3), respectively, and solid lines represent the exact solutions. The used parameters are 
S7 = Sg = 3 X 10^, Si3 = 10^, other collision parameters are 10^. The relative errors of 
the density, pressure, velocity and temperature for case 1 are 0.398%, 0.205%, 0.592% and 
0.310%, respectively. The relative errors for case 2 are 0.344%, 0.130%, 0.408% and 0.287%, 
respectively. And the relative errors for case 3 are 0.334%, 0.117%, 0.372% and 0.283%, 
respectively. 

(c) Colella explosion wave 
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FIG. 4: The MRT LB numerical and exact solutions for Lax shock tube at time t = 0.2. 



(23) 



For the problem, the initial condition is 

{p,ui,U2,T)\l = (LO, 0.0, 0.0, 1000.0), x < 0. 
{p,UuU2,T)\ii = (LO, 0.0, 0.0, 0.01), x > 0. 

This is a strong temperature discontinuity problem that can be used to study the ro- 
bustness and precision of numerical methods. Figure 5 gives density, pressure, velocity and 
temperature results at t = 0.1. The red squares and green circles correspond to simulation 
results with different grid sizes and time steps: dx = dy = 0.002, dt = 2 x 10~^ (case 1), 
and dx = dy = 0.001, dt = 10^® (case 2), respectively, and solid lines represent the exact 
solutions. Here, the parameters are sy = ss = 5 x 10"^, sn = sis = 5 x 10^, other values 
of s adopt 10^. The relative errors of the density, pressure, velocity and temperature for 
case 1 are 1.69%, 1.11%, 1.60% and 0.779%, respectively. The relative errors for case 2 are 
1.68%, 1.11%, 1.59% and 0.777%, respectively. The oscillations at the interface are difficult 
to eliminate completely in our simulations. 

(d) Collision of two strong shocks 

The initial condition can be write as: 

{p,Ui,U2,T)\l = (5.99924,19.5975,0.0,76.8254), x < 0. 
{p,Ui,U2,T)\r = (5.99242,-6.19633,0.0,7.69222), x > 0. 



(24) 
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FIG. 5: The MRT simulation results and exact solutions for the Colella explosion wave at time 
t = 0.1. 

The MRT and SRT numerical results and exact solutions at time t = 0.1 are shown in 
Fig.6, where the common parameters are dx = dy = 0.003, dt = 10~^, the collision matrix 
in MRT is S5 = Sg = 5 X 10'^, S7 = Sg = 3 x lO'*, other values of s are 10^, and the relaxation 
time in SRT is r = 10~^. The red squares and green circles correspond to the MRT and 
SRT simulation results, respectively, and solid lines represent the exact solutions. Compared 
with the simulation results of SRT, we can find that the oscillations at the discontinuity are 
weaker in MRT simulation. 

(e) High Mach number shock 

In order to test the stability of the new model, we construct a new shock tube problem 
with the initial condition, 

(p,Mi,M2,T)|i = (5.0,45.0,0.0,10.0), x < 0. 

[25) 

{p,Ui,U2,T)\r = (6.0, -20.0,0.0,5.0), x > 0. 

The Mach number of the left side is 10.1, and the right is 6.3. And this test is successfully 
passed by the MRT LB, but failed by the SRT. Figure 7 shows comparison of the MRT 
LB results and exact solutions at t = 0.018, where the parameters are dx = dy = 0.003, 
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FIG. 6: The MRT and SRT numerical results and exact solution for collision of two strong shocks 
at time t = 0.1. 

dt = 10^^, S5 = sg = 1.5 X 10^, Sio = 5 X 10^, other values of s are 10^. Circle symbols 
correspond to MRT simulation results, and solid lines represent the exact solutions. 

B. Shock reflections 

We adopt the macroscopic variable boundary conditions in Figs. 8-11. The values of 
the distribution functions on boundaries are assigned with the corresponding values of the 
equilibrium distribution functions. The determination methods of macroscopic quantities 
are explained combining with specific problems. Here we study two gas dynamic problems: 
regular and Mach shock reflections. 

(a) Regular shock reflection 

In the problem, we have performed a 30° shock reflection, and the corresponding Mach 
number is 5. The computational domain is a rectangle with length 0.6 and height 0.2, which 
is divided into 300 x 100. The reflecting wall lies at the bottom of the domain (reflecting 
boundary condition denotes the y component of the fluid velocity on the boundary is reverse 
to that of interior point), and the linear extrapolation technique is applied to define the 
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FIG. 7: MRT LB results and exact solutions for the high Mach number shock tube problem at 
t = 0.018. 

values of the macroscopic quantities on the right-hand boundary. The other two sides adopt 
the Dirichlet boundary conditions: 

(p, «i, U2, T)\o,y,t = (1.0, 5.0, 0.0, 0.5), ^^^^ 
(p, «i, M2, T)|,,o.2,t = (2.27273, 4.3, -1.21244, 1.76). 

Initially, the entire interior zone is set the values of the left boundary. Parameters in the 
simulation are as follows: dx = dy = 0.002, dt = 10~^, S5 = 10^, se 2 x 10^, S7 — Ss — 10^, 
other collision parameters arc 10^. Figure 8 shows the pressure contour at time t = 0.3, and 
the density, temperature contours have similar results. From black to white, the grey level 
corresponds to the increase of the pressure. The result shows that the new MRT model has 
the ability to accurately capture the shock front, 
(b) Mach reflection problem 

This problem is on an unsteady shock reflection. A planar shock impacts an oblique 
surface which is at a 30° angle to the propagation direction of the shock. The fluid in front 
of the shock has zero velocity, and the shock Mach number is 1.5. The initial condition is 
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FIG. 8: Pressure contour of regular shock reflection at time t = 0.3. 



as follows: 



{p,Ui,U2,T) 



x,y,0~ 



(3.17647, 0.555556 cos30°, -0.555556 sin 30°, 0.839506), ii y > h{x,0) 

(2.0,0.0,0.0,0.5), iiy<h{x,0) 

(27) 

where h{x, t) = tan60°(a; — 150dx) — 1.5t/ sin 30°. The computational domain is divided into 
600 X 300. At the bottom boundary, reflecting boundary condition is used from the 150th 
grid, and the left side adopts the values of the initial post-shock flow; the left boundary is also 
assigned values of the initial post-shock flow, and at the right boundary the extrapolation 
technique is applied; at the top boundary, the macroscopic variables are assigned using the 
same method of the right boundary when x > g{t), and are set the same values as the left 
boundary's when x < g(t), where g{t) = 150dx + tan30°(0.9 + 1.5t/ sin 30°), dx is the grid 
size in simulation. In Fig. 9 we show the results of density, pressure, velocity in x direction 
and temperature contours at t = 0.25 in the part of [50,450] x [0,200]. Parameters in the 



simulation are dx = dy = 0.003, dt 



10- 



S5 



10^ S6 = 5 X 102, S7 



S8 



10^, and 



other collision parameters are 10^. The simulation results are accordant with those of other 



numerical methods 



471, 



51 



-l53|. 



C. Shock wave reaction on cylindrical bubble problem 

The problems are as follows: A planar shock wave with the Mach number 1.22 impinges 
on a cylindrical bubble with different densities. In the first case the bubble has a lower 
density. In the second case the bubble's density is higher. 
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FIG. 9: Mach reflection of shock wave. Figs, (a)-(d) shows the density, pressure, velocity in x 
direction, and temperature at time t = 0.25, respectively. From black to white, the grey level 
corresponds to the increase of values. 

Initial conditions for the first case are 



{p,Ui,U2,p) U,2/,0= < 



and for the second case are 



{p,Ui,U2,p) 



x,y,0- 



(1, 0, 0, 1) , pre — shock, 

(1.28, -0.3774, 0, 1.6512) , post - shock, 
(0.1358, 0, 0, 1) , bubble, 

(1, 0, 0, 1) , pre — shock, 

(1.28, -0.3774, 0, 1.6512) , post - shock, 
(3.1538,0,0,1), bubble. 



(28) 



(29) 



In the simulations, the right side adopts the values of the initial post-shock fiow; the 
extrapolation technique is applied at the left boundary, and refiecting boundary conditions 
are imposed on the top and bottom. The common parameters arc as follows: dx = dy = 
0.003, dt = 10~^. When simulate the low density cylindrical bubble, the collision parameters 
are — se — — ss — 10^, and others are 10^; when simulate the high density bubble, 
the coUision parameters are S5 — 10^, and s — 10^ for the others. In Fig. 10(a), from top to 
bottom, the three plots show the density contours at the times t — 0, 0.5, 0.65, respectively. 
In Fig. 10(b), from top to bottom, the three plots show the density contours at the times 



20 



50 



50 100 150 200 250 300 50 100 150 200 250 300 




250 300 



FIG. 10: Snapshots of shock wave reaction on single bubble. The left column is for the process 
with initial condition (I28p , and the right column is for the process with initial condition ()29|) . From 
black to white, the density value increases. 

t = 0, 0.6, 0.9, respectively. These results are accordant with those from other methods 
and experiment |55|] . The surface of bubbles is comparatively smooth, which indicates that 
the MRT modle has high accuracy and resolution. 





D. Couette flow 



In order to demonstrate the accuracy of the model, numerical simulations of the incom- 
pressible Couette flow are carried out. Consider a viscous fluid flow between two infinite 
parallel flat plates, separated by a distance of D. The initial state of the fluid is p = 1, 
T = 1, U = 0. At time t = the two plates start to move at velocities U, —U, respectively, 
{U = 0.1). The periodic boundary condition is adopted in the x direction. The top and bot- 
tom boundaries are constant speed and constant temperature boundaries {U = 0.1, T = 1). 
The analytical solution of horizontal velocity along a vertical line is as follows: 

u = 2yU/D - — exp[ -^t] sm{—y), 

where j is an integer, the two walls locate at y = ±D/2. 

We carried out a set of simulations: dx = dy = 0.004, dt = 10~^, NX x NY = 16 x 32 
(case 1), dx = dy = 0.002, dt = 5 x 10-^ NX x NY = 32 x 64(case 2), dx = dy = 0.001, 
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FIG. 11: Comparison of the MRT LB simulation results and exact solutions for the horizontal 
velocity distribution of Couette how at time t = 57.5. 

dt = 2.5 X 10"^ NX X NY = 64 x 128 (case 3). All of the collision parameters are 10^ 
Figure 11 shows a comparison of the MRT LB simulation results and exact solution for the 
horizontal velocity distribution at time t = 57.5. The black squares, red triangles and green 
circles correspond to simulation results of case 1, case 2, and case 3, respectively, and solid 
line represents the exact solution. The relative errors of the horizontal velocity for the three 
cases are 7.02%, 3.86% and 2.06%, respectively. So this model is of first order accuracy. 



VI. CONCLUSION 



An energy- conserving multiple-relaxation-time finite difference lattice Boltzmann model 
for compressible flows is proposed. This model is based on a 16-discrete-velocity model de- 
signed by Kataoka and Tsutaharaj4l|. The collision step is flrst calculated in the moment 
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space and then mapped back to the velocity space. The moment space and correspond- 
ing transformation matrix are constructed according to the group representation theory. 
Equihbria of the nonconserved moments are chosen according to the need of recovering 
compressible Navier-Stokes equations through the Chapman-Enskog expansion. In the new 
model different transport coefficients, such as viscosity and heat conductivity, are related 
to different collision parameters. Consequently, they can be controlled independently. This 
flexibility makes the MRT model contain more physical information and easier to satisfy the 
von Neumann stability condition than its SRT counterpart. The new model passed well- 
known benchmark tests, including (i) shock tubes, such as the Sod, Lax, Colella explosion 
wave, collision of two strong shocks and a new shock tube with high Mach number, (ii) regu- 
lar and Mach shock reflections, (iii) shock wave reaction on cylindrical bubble problems, and 
(iv) Couette flow. This model works for both low and high speeds compressible flows. Both 
the LB and traditional CFD approach work when the Knudsen number is very small, but in 
the vicinity of shock wave, the system is in a nonequilibrium state, and the traditional Euler 
and Navier-Stokes descriptions are problematic. For such problems, LB has more sound 
physical ground. 
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Appendix. Spatial and temporal discretization effects 

In our simulations, the time evolution is based on the usual first-order upwind scheme, 
while space discretization is performed through a Lax-Wendroff scheme. 

^Ji-m + St)-Mt)]/St, (30a) 
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'^i-yfi = Vi^[fi{x + 5x,y) - fi{x - 5x,y)]/2Sx 
+Viy[fiix, y + 6y)- fi{x, y - 6y)]/25y 
-VixSt[fi{x + 6x, y) - 2fi{x, y) + f^{x - 6x, y)]/26x^ 
-vl6t[fi{x,y + 6y) - 2f,{x,y) + f,{x,y - 6y)]/26y'. (30b) 

If we perform the following series expansion up to second order in Eq.f l30ap and Eq.f lSObp 

Mt + 5t) = m + + l^t'^^f,, (31a) 

d 1 (9^ 

fi{x + 5x, y, t) = fi{x, y, t) + Sx—fi + l^^^'^^ft^ (31b) 

d I d"^ 

fi{x - 5x, y, t) = fi{x, y, t) - Sx—fi + 2^^^g^f'' (^-^^) 

we can get 

d (9 1 (9^ 

V. . V/. ^ + - '^^n - "J^p., (32b) 

So a more accurate LB equation solved using the updating schemes mentioned above is 



To investigate the effect of the supplementary terms in Eq.( !33l) . a similar Chapman- 
Enskog expansion procedure may be considered. The zeroth and first order LB equations 
Eq.f l7a|) . Eq. (l7bp remain unchanged when performing the Chapman-Enskog expansion, while 
the second-order equation (ITcI) becomes, 

^ ^(0) , d (1) 6t d , ^, d d (i) 

It can be converted into moment space to obtain: 

9 m, d ^(1) 6t d ^ ^(l)^ , d d ^ ^(1) 



^ d_ ^ ?(0) _ ^ J.2 

where = M(fQ,I)M~^. From the equation, we obtain 



F F - ^(0) "^p.2_^?{0)_ f(2) ("jrN 

2^'''^''d^^d^J^ ~^^^-dxf^^ ' ^^^^ 



^ + 5t7^7^/r = 0, (36a) 
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In this way the recovered NS equations follows: 

dp , <9j^. djy r+ ^ "9 

(37b) 
(37c) 

| + |l(. + PW./p] + |l(e + Pb-» 

d , dT duy du^ duy du^ 

= — Ai(2— + Uy—- + u^— Ux^r- + Uy^—)\ 

ox ox ox ox oy oy 

d , dT du^ du^ duy duy 

-^2(2^ + ^y^r- + + a ) 

oy oy oy ox ox oy 

d d 1 

where 

p's = pRTi- + f ), p: = pRT{- + f ), a; = pR^n^ + 1), a; = pR^n^ + 1). 

S5 2 Se 2 S7 2 ss 2 

When 6t approaches 0, equations fl37ap - (]37dp reduce to the Eqs. fll3a|) - (]13dp . 
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